function Rates = my_compute_rate(qbit,Nsub,L,TRX,fc,BW_des,Ptx,showfig)
% qbit=1;

% Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
% L = 15e-1;
% 
% factor = 0.3;
% x1 = 40e-1*factor;
% y1 = 0e-1*factor;
% z1 = 90e-1*factor;
% x2 = -18e-1*factor;
% y2 = 0e-1*factor;
% z2 = 20e-1*factor;
x1 = TRX(1);y1 = TRX(2);z1 = TRX(3);
x2 = TRX(4);y2 = TRX(5);z2 = TRX(6);
z12 = z1^2;z22 = z2^2;
% 
% fc = 30e9;
c = 3e8;
% BW_des = 1.5e9;
% Ptx = 0;%dBm


Ptx_lin = 10.^((Ptx-30)/10); % in watts
Lb = c/fc;
spac = Lb/2;
hp = 6.625e-34;
kb = 1.38e-23;
Tk = 290;


LIS_f = 0;
N = round(L/spac);

xcoord=linspace(-L/2,L/2,N);
ycoord=xcoord;
[xx,yy]=meshgrid(xcoord,ycoord);
a = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2))/2;
a_vec = linspace(min(a(:)),max(a(:)),Ngrid)';% column vector
Distmat = a*2;
distgrid = a_vec*2;
Distmat1 = sqrt((xx-x1).^2+(yy-y1).^2+z1^2);
Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);

[ModBFmat,QBFmat,modfn] = Beamforming(qbit,L,TRX,fc,BW_des,Ptx,showfig);


        freqvec=2*pi*linspace(fc-BW_des/2,fc+BW_des/2,Nsub);
        antt_delays=Distmat(:)/c;
        AF_product=antt_delays*freqvec;
        h_0 = (Distmat1(:)*freqvec/c).*(Distmat2(:)*freqvec/c);
        
        Wb_mch=exp(-1i*AF_product)./h_0;
        SCpowers_new_1=abs(ModBFmat(:).'*Wb_mch).^2;
        SCpowers_old_1=abs(QBFmat(:).'*Wb_mch).^2;
        SCpowers_ideal = zeros(size(freqvec));
        SCpowers_ideal(abs(freqvec/2/pi-fc)<=BW_des/2) = 1;
        Pow = sum(modfn.^2)*spac^-4*(a_vec(end)-a_vec(1))/Ngrid*c/(freqvec(2)-freqvec(1))*pi./h_0(end/2,end/2)^2;
        temp = [Pow sum(SCpowers_old_1)];
        SCpowers_ideal_1 = SCpowers_ideal*Pow/sum(SCpowers_ideal);
       % SCpowers_ideal_1 = max(mean(SCpowers_old_1),mean(SCpowers_new_1))*ones(size(SCpowers_new));
 
%         figure;
%         plot(freqvec/2/pi,10*log10([SCpowers_old' SCpowers_new']))
        %% Compute rate
        
        Delta_f=freqvec(2)-freqvec(1);
        Rates = zeros(2,1);
        phin_nr=hp*freqvec; % ε=hν
        phin_dr=exp(hp*freqvec/(kb*Tk))-1;% Bose distribution factor
        
        phin=phin_nr./phin_dr;  % in watts
        
        Hwb_sq=SCpowers_old_1;%
        N_eff=phin./Hwb_sq;
        % mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        % powalloc=mu_df-N_eff;
        % powalloc(powalloc<0)=0;
        powalloc = Ptx_lin/BW_des;
        SNR_eff=powalloc./N_eff;
        Rates(1)=Delta_f*sum(log2(1+SNR_eff));
        
        Hwb_sq=SCpowers_new_1;
        N_eff=phin./Hwb_sq;
        % mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        % powalloc=mu_df-N_eff;
        % powalloc(powalloc<0)=0;
        powalloc = Ptx_lin/BW_des;
        SNR_eff=powalloc./N_eff;
        Rates(2)=Delta_f*sum(log2(1+SNR_eff));
        
        Hwb_sq=SCpowers_ideal_1;%
        N_eff=phin./Hwb_sq;
        % mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        % powalloc=mu_df-N_eff;
        % powalloc(powalloc<0)=0;
        powalloc = Ptx_lin/BW_des;
        SNR_eff=powalloc./N_eff;
        Rates(3)=Delta_f*sum(log2(1+SNR_eff));

        
        
        